Cardiac mechanics Benchmark - Problem 2

Cardiac mechanics Benchmark - Problem 2#

This code implements problem 2 in the Cardiac Mechanic Benchmark paper

Land S, Gurev V, Arens S, Augustin CM, Baron L, Blake R, Bradley C, Castro S, Crozier A, Favino M, Fastl TE. Verification of cardiac mechanics software: benchmark problems and solutions for testing active and passive material behaviour. Proc. R. Soc. A. 2015 Dec 8;471(2184):20150641.

from pathlib import Path
import dolfin
try:
    from dolfin_adjoint import Constant, DirichletBC, Mesh, interpolate
except ImportError:
    from dolfin import DirichletBC, Constant, interpolate, Mesh
import pulse
import cardiac_geometries as cg
from fenics_plotly import plot
geo_path = Path("geometry")
if not geo_path.is_dir():
    cg.create_benchmark_geometry_land15(outdir=geo_path)
geo = cg.geometry.Geometry.from_folder(geo_path)
geometry = pulse.HeartGeometry(
    mesh=geo.mesh,
    markers=geo.markers,
    marker_functions=pulse.MarkerFunctions(vfun=geo.vfun, ffun=geo.ffun),
    microstructure=pulse.Microstructure(f0=geo.f0, s0=geo.s0, n0=geo.n0),
)
Info    : Meshing 1D...
Info    : [  0%] Meshing curve 1 (Line)
Info    : [ 10%] Meshing curve 2 (Line)
Info    : [ 10%] Meshing curve 3 (Ellipse)
Info    : [ 20%] Meshing curve 4 (Ellipse)
Info    : [ 20%] Meshing curve 7 (Ellipse)
Info    : [ 30%] Meshing curve 8 (Line)
Info    : [ 30%] Meshing curve 9 (Ellipse)
Info    : [ 40%] Meshing curve 12 (Circle)
Info    : [ 40%] Meshing curve 16 (Circle)
Info    : [ 50%] Meshing curve 24 (Ellipse)
Info    : [ 50%] Meshing curve 25 (Line)
Info    : [ 60%] Meshing curve 26 (Ellipse)
Info    : [ 60%] Meshing curve 29 (Circle)
Info    : [ 70%] Meshing curve 33 (Circle)
Info    : [ 70%] Meshing curve 41 (Ellipse)
Info    : [ 80%] Meshing curve 42 (Line)
Info    : [ 80%] Meshing curve 43 (Ellipse)
Info    : [ 90%] Meshing curve 46 (Circle)
Info    : [ 90%] Meshing curve 50 (Circle)
Info    : [100%] Meshing curve 63 (Circle)
Info    : [100%] Meshing curve 67 (Circle)
Info    : Done meshing 1D (Wall 0.0197534s, CPU 0.016493s)
Info    : Meshing 2D...
Info    : [  0%] Meshing surface 1 (Plane, Frontal-Delaunay)
Info    : [ 10%] Meshing surface 13 (Surface, Frontal-Delaunay)
Info    : [ 20%] Meshing surface 17 (Surface, Frontal-Delaunay)
Info    : [ 20%] Meshing surface 20 (Surface, Frontal-Delaunay)
Info    : [ 30%] Meshing surface 21 (Plane, Frontal-Delaunay)
Info    : [ 40%] Meshing surface 30 (Surface, Frontal-Delaunay)
Info    : [ 40%] Meshing surface 34 (Surface, Frontal-Delaunay)
Info    : [ 50%] Meshing surface 37 (Surface, Frontal-Delaunay)
Info    : [ 50%] Meshing surface 38 (Plane, Frontal-Delaunay)
Info    : [ 60%] Meshing surface 47 (Surface, Frontal-Delaunay)
Info    : [ 70%] Meshing surface 51 (Surface, Frontal-Delaunay)
Info    : [ 70%] Meshing surface 54 (Surface, Frontal-Delaunay)
Info    : [ 80%] Meshing surface 55 (Plane, Frontal-Delaunay)
Info    : [ 90%] Meshing surface 64 (Surface, Frontal-Delaunay)
Info    : [ 90%] Meshing surface 68 (Surface, Frontal-Delaunay)
Info    : [100%] Meshing surface 71 (Surface, Frontal-Delaunay)
Info    : Done meshing 2D (Wall 0.0214609s, CPU 0.021555s)
Info    : Meshing 3D...
Info    : 3D Meshing 4 volumes with 1 connected component
Info    : Tetrahedrizing 599 nodes...
Info    : Done tetrahedrizing 607 nodes (Wall 0.0058539s, CPU 0.005889s)
Info    : Reconstructing mesh...
Info    :  - Creating surface mesh
Info    :  - Identifying boundary edges
Info    :  - Recovering boundary
Info    : Done reconstructing mesh (Wall 0.0127351s, CPU 0.009543s)
Info    : Found volume 2
Info    : Found volume 1
Info    : Found volume 4
Info    : Found volume 3
Info    : It. 0 - 0 nodes created - worst tet radius 0.904667 (nodes removed 0 0)
Info    : 3D refinement terminated (599 nodes total):
Info    :  - 0 Delaunay cavities modified for star shapeness
Info    :  - 0 nodes could not be inserted
Info    :  - 1934 tetrahedra created in 4.1117e-05 sec. (47036505 tets/s)
Info    : Done meshing 3D (Wall 0.0213887s, CPU 0.018004s)
Info    : Optimizing mesh...
Info    : Optimizing volume 1
Info    : Optimization starts (volume = 805.7) with worst = 0.0233132 / average = 0.739813:
Info    : 0.00 < quality < 0.10 :         7 elements
Info    : 0.10 < quality < 0.20 :         5 elements
Info    : 0.20 < quality < 0.30 :         5 elements
Info    : 0.30 < quality < 0.40 :         7 elements
Info    : 0.40 < quality < 0.50 :         7 elements
Info    : 0.50 < quality < 0.60 :        45 elements
Info    : 0.60 < quality < 0.70 :        59 elements
Info    : 0.70 < quality < 0.80 :       127 elements
Info    : 0.80 < quality < 0.90 :       188 elements
Info    : 0.90 < quality < 1.00 :        36 elements
Info    : 16 edge swaps, 0 node relocations (volume = 805.7): worst = 0.302779 / average = 0.760591 (Wall 0.000340498s, CPU 0.000383s)
Info    : No ill-shaped tets in the mesh :-)
Info    : 0.00 < quality < 0.10 :         0 elements
Info    : 0.10 < quality < 0.20 :         0 elements
Info    : 0.20 < quality < 0.30 :         0 elements
Info    : 0.30 < quality < 0.40 :         8 elements
Info    : 0.40 < quality < 0.50 :         8 elements
Info    : 0.50 < quality < 0.60 :        45 elements
Info    : 0.60 < quality < 0.70 :        59 elements
Info    : 0.70 < quality < 0.80 :       128 elements
Info    : 0.80 < quality < 0.90 :       188 elements
Info    : 0.90 < quality < 1.00 :        36 elements
Info    : Optimizing volume 2
Info    : Optimization starts (volume = 805.7) with worst = 0.0233132 / average = 0.736103:
Info    : 0.00 < quality < 0.10 :         7 elements
Info    : 0.10 < quality < 0.20 :         5 elements
Info    : 0.20 < quality < 0.30 :         5 elements
Info    : 0.30 < quality < 0.40 :         8 elements
Info    : 0.40 < quality < 0.50 :        11 elements
Info    : 0.50 < quality < 0.60 :        44 elements
Info    : 0.60 < quality < 0.70 :        59 elements
Info    : 0.70 < quality < 0.80 :       125 elements
Info    : 0.80 < quality < 0.90 :       190 elements
Info    : 0.90 < quality < 1.00 :        33 elements
Info    : 16 edge swaps, 0 node relocations (volume = 805.7): worst = 0.302779 / average = 0.754569 (Wall 0.000297568s, CPU 0.000308s)
Info    : No ill-shaped tets in the mesh :-)
Info    : 0.00 < quality < 0.10 :         0 elements
Info    : 0.10 < quality < 0.20 :         0 elements
Info    : 0.20 < quality < 0.30 :         0 elements
Info    : 0.30 < quality < 0.40 :         8 elements
Info    : 0.40 < quality < 0.50 :        14 elements
Info    : 0.50 < quality < 0.60 :        45 elements
Info    : 0.60 < quality < 0.70 :        61 elements
Info    : 0.70 < quality < 0.80 :       126 elements
Info    : 0.80 < quality < 0.90 :       187 elements
Info    : 0.90 < quality < 1.00 :        33 elements
Info    : Optimizing volume 3
Info    : Optimization starts (volume = 805.689) with worst = 0.0309823 / average = 0.745513:
Info    : 0.00 < quality < 0.10 :         7 elements
Info    : 0.10 < quality < 0.20 :         3 elements
Info    : 0.20 < quality < 0.30 :         2 elements
Info    : 0.30 < quality < 0.40 :         6 elements
Info    : 0.40 < quality < 0.50 :        11 elements
Info    : 0.50 < quality < 0.60 :        42 elements
Info    : 0.60 < quality < 0.70 :        61 elements
Info    : 0.70 < quality < 0.80 :       128 elements
Info    : 0.80 < quality < 0.90 :       188 elements
Info    : 0.90 < quality < 1.00 :        33 elements
Info    : 12 edge swaps, 0 node relocations (volume = 805.689): worst = 0.319724 / average = 0.76101 (Wall 0.000367529s, CPU 0.000341s)
Info    : No ill-shaped tets in the mesh :-)
Info    : 0.00 < quality < 0.10 :         0 elements
Info    : 0.10 < quality < 0.20 :         0 elements
Info    : 0.20 < quality < 0.30 :         0 elements
Info    : 0.30 < quality < 0.40 :         6 elements
Info    : 0.40 < quality < 0.50 :        11 elements
Info    : 0.50 < quality < 0.60 :        42 elements
Info    : 0.60 < quality < 0.70 :        61 elements
Info    : 0.70 < quality < 0.80 :       132 elements
Info    : 0.80 < quality < 0.90 :       187 elements
Info    : 0.90 < quality < 1.00 :        32 elements
Info    : Optimizing volume 4
Info    : Optimization starts (volume = 805.685) with worst = 0.0552945 / average = 0.747184:
Info    : 0.00 < quality < 0.10 :         5 elements
Info    : 0.10 < quality < 0.20 :         5 elements
Info    : 0.20 < quality < 0.30 :         1 elements
Info    : 0.30 < quality < 0.40 :         7 elements
Info    : 0.40 < quality < 0.50 :        12 elements
Info    : 0.50 < quality < 0.60 :        35 elements
Info    : 0.60 < quality < 0.70 :        66 elements
Info    : 0.70 < quality < 0.80 :       129 elements
Info    : 0.80 < quality < 0.90 :       186 elements
Info    : 0.90 < quality < 1.00 :        34 elements
Info    : 11 edge swaps, 0 node relocations (volume = 805.685): worst = 0.326537 / average = 0.763665 (Wall 0.000487344s, CPU 0s)
Info    : No ill-shaped tets in the mesh :-)
Info    : 0.00 < quality < 0.10 :         0 elements
Info    : 0.10 < quality < 0.20 :         0 elements
Info    : 0.20 < quality < 0.30 :         0 elements
Info    : 0.30 < quality < 0.40 :         7 elements
Info    : 0.40 < quality < 0.50 :        11 elements
Info    : 0.50 < quality < 0.60 :        34 elements
Info    : 0.60 < quality < 0.70 :        64 elements
Info    : 0.70 < quality < 0.80 :       133 elements
Info    : 0.80 < quality < 0.90 :       187 elements
Info    : 0.90 < quality < 1.00 :        33 elements
Info    : Done optimizing mesh (Wall 0.00424632s, CPU 0.002473s)
Info    : Optimizing mesh (Netgen)...
Info    : Optimizing volume 1
Info    : CalcLocalH: 187 Points 472 Elements 370 Surface Elements 
Info    : Remove Illegal Elements 
Info    : 169 illegal tets 
Info    : SplitImprove 
Info    : badmax = 21.389 
Info    : 16 splits performed 
Info    : SwapImprove  
Info    : 16 swaps performed 
Info    : SwapImprove2  
Info    : 1 swaps performed 
Info    : 124 illegal tets 
Info    : SplitImprove 
Info    : badmax = 99.5365 
Info    : 17 splits performed 
Info    : SwapImprove  
Info    : 7 swaps performed 
Info    : SwapImprove2  
Info    : 0 swaps performed 
Info    : 90 illegal tets 
Info    : SplitImprove 
Info    : badmax = 99.5365 
Info    : 15 splits performed 
Info    : SwapImprove  
Info    : 4 swaps performed 
Info    : SwapImprove2  
Info    : 0 swaps performed 
Info    : 64 illegal tets 
Info    : SplitImprove 
Info    : badmax = 99.5365 
Info    : 15 splits performed 
Info    : SwapImprove  
Info    : 5 swaps performed 
Info    : SwapImprove2  
Info    : 0 swaps performed 
Info    : 29 illegal tets 
Info    : SplitImprove 
Info    : badmax = 101.081 
Info    : 8 splits performed 
Info    : SwapImprove  
Info    : 0 swaps performed 
Info    : SwapImprove2  
Info    : 0 swaps performed 
Info    : 6 illegal tets 
Info    : SplitImprove 
Info    : badmax = 101.081 
Info    : 2 splits performed 
Info    : SwapImprove  
Info    : 0 swaps performed 
Info    : SwapImprove2  
Info    : 1 swaps performed 
Info    : 0 illegal tets 
Info    : Volume Optimization 
Info    : CombineImprove 
Info    : 19 elements combined 
Info    : ImproveMesh 
Info    : Total badness = 1449.12 
Info    : Total badness = 1358.97 
Info    : SplitImprove 
Info    : badmax = 58.0222 
Info    : 2 splits performed 
Info    : ImproveMesh 
Info    : Total badness = 1367.36 
Info    : Total badness = 1360.57 
Info    : SwapImprove  
Info    : 65 swaps performed 
Info    : SwapImprove2  
Info    : 0 swaps performed 
Info    : ImproveMesh 
Info    : Total badness = 1229.14 
Info    : Total badness = 1167.69 
Info    : CombineImprove 
Info    : 10 elements combined 
Info    : ImproveMesh 
Info    : Total badness = 1070.2 
Info    : Total badness = 1065.73 
Info    : SplitImprove 
Info    : badmax = 15.4382 
Info    : 0 splits performed 
Info    : ImproveMesh 
Info    : Total badness = 1065.73 
Info    : Total badness = 1065.54 
Info    : SwapImprove  
Info    : 14 swaps performed 
Info    : SwapImprove2  
Info    : 0 swaps performed 
Info    : ImproveMesh 
Info    : Total badness = 1051.45 
Info    : Total badness = 1039.94 
Info    : CombineImprove 
Info    : 2 elements combined 
Info    : ImproveMesh 
Info    : Total badness = 1021.92 
Info    : Total badness = 1021.4 
Info    : SplitImprove 
Info    : badmax = 15.4382 
Info    : 0 splits performed 
Info    : ImproveMesh 
Info    : Total badness = 1021.4 
Info    : Total badness = 1021.37 
Info    : SwapImprove  
Info    : 7 swaps performed 
Info    : SwapImprove2  
Info    : 0 swaps performed 
Info    : ImproveMesh 
Info    : Total badness = 1018.11 
Info    : Total badness = 1013.1 
Info    : Optimizing volume 2
Info    : CalcLocalH: 187 Points 474 Elements 370 Surface Elements 
Info    : Remove Illegal Elements 
Info    : 170 illegal tets 
Info    : SplitImprove 
Info    : badmax = 21.389 
Info    : 17 splits performed 
Info    : SwapImprove  
Info    : 22 swaps performed 
Info    : SwapImprove2  
Info    : 1 swaps performed 
Info    : 119 illegal tets 
Info    : SplitImprove 
Info    : badmax = 99.5365 
Info    : 15 splits performed 
Info    : SwapImprove  
Info    : 5 swaps performed 
Info    : SwapImprove2  
Info    : 0 swaps performed 
Info    : 85 illegal tets 
Info    : SplitImprove 
Info    : badmax = 99.5365 
Info    : 14 splits performed 
Info    : SwapImprove  
Info    : 5 swaps performed 
Info    : SwapImprove2  
Info    : 1 swaps performed 
Info    : 54 illegal tets 
Info    : SplitImprove 
Info    : badmax = 433.677 
Info    : 14 splits performed 
Info    : SwapImprove  
Info    : 2 swaps performed 
Info    : SwapImprove2  
Info    : 1 swaps performed 
Info    : 27 illegal tets 
Info    : SplitImprove 
Info    : badmax = 433.677 
Info    : 9 splits performed 
Info    : SwapImprove  
Info    : 1 swaps performed 
Info    : SwapImprove2  
Info    : 0 swaps performed 
Info    : 0 illegal tets 
Info    : Volume Optimization 
Info    : CombineImprove 
Info    : 16 elements combined 
Info    : ImproveMesh 
Info    : Total badness = 1477.7 
Info    : Total badness = 1376.18 
Info    : SplitImprove 
Info    : badmax = 58.02 
Info    : 2 splits performed 
Info    : ImproveMesh 
Info    : Total badness = 1384.57 
Info    : Total badness = 1378.78 
Info    : SwapImprove  
Info    : 62 swaps performed 
Info    : SwapImprove2  
Info    : 0 swaps performed 
Info    : ImproveMesh 
Info    : Total badness = 1270.02 
Info    : Total badness = 1206.98 
Info    : CombineImprove 
Info    : 11 elements combined 
Info    : ImproveMesh 
Info    : Total badness = 1099.72 
Info    : Total badness = 1091.5 
Info    : SplitImprove 
Info    : badmax = 17.0774 
Info    : 0 splits performed 
Info    : ImproveMesh 
Info    : Total badness = 1091.5 
Info    : Total badness = 1090.85 
Info    : SwapImprove  
Info    : 22 swaps performed 
Info    : SwapImprove2  
Info    : 0 swaps performed 
Info    : ImproveMesh 
Info    : Total badness = 1068.29 
Info    : Total badness = 1043.93 
Info    : CombineImprove 
Info    : 3 elements combined 
Info    : ImproveMesh 
Info    : Total badness = 1020.02 
Info    : Total badness = 1018.36 
Info    : SplitImprove 
Info    : badmax = 15.4382 
Info    : 0 splits performed 
Info    : ImproveMesh 
Info    : Total badness = 1018.36 
Info    : Total badness = 1018.28 
Info    : SwapImprove  
Info    : 7 swaps performed 
Info    : SwapImprove2  
Info    : 0 swaps performed 
Info    : ImproveMesh 
Info    : Total badness = 1015.73 
Info    : Total badness = 1012.33 
Info    : Optimizing volume 3
Info    : CalcLocalH: 187 Points 471 Elements 370 Surface Elements 
Info    : Remove Illegal Elements 
Info    : 170 illegal tets 
Info    : SplitImprove 
Info    : badmax = 14.227 
Info    : 18 splits performed 
Info    : SwapImprove  
Info    : 21 swaps performed 
Info    : SwapImprove2  
Info    : 0 swaps performed 
Info    : 121 illegal tets 
Info    : SplitImprove 
Info    : badmax = 35.9306 
Info    : 16 splits performed 
Info    : SwapImprove  
Info    : 5 swaps performed 
Info    : SwapImprove2  
Info    : 0 swaps performed 
Info    : 88 illegal tets 
Info    : SplitImprove 
Info    : badmax = 35.9306 
Info    : 17 splits performed 
Info    : SwapImprove  
Info    : 6 swaps performed 
Info    : SwapImprove2  
Info    : 1 swaps performed 
Info    : 49 illegal tets 
Info    : SplitImprove 
Info    : badmax = 87.3503 
Info    : 13 splits performed 
Info    : SwapImprove  
Info    : 3 swaps performed 
Info    : SwapImprove2  
Info    : 0 swaps performed 
Info    : 16 illegal tets 
Info    : SplitImprove 
Info    : badmax = 87.3503 
Info    : 5 splits performed 
Info    : SwapImprove  
Info    : 1 swaps performed 
Info    : SwapImprove2  
Info    : 0 swaps performed 
Info    : 3 illegal tets 
Info    : SplitImprove 
Info    : badmax = 87.3503 
Info    : 1 splits performed 
Info    : SwapImprove  
Info    : 0 swaps performed 
Info    : SwapImprove2  
Info    : 0 swaps performed 
Info    : 0 illegal tets 
Info    : Volume Optimization 
Info    : CombineImprove 
Info    : 23 elements combined 
Info    : ImproveMesh 
Info    : Total badness = 1314.97 
Info    : Total badness = 1256.73 
Info    : SplitImprove 
Info    : badmax = 57.7433 
Info    : 1 splits performed 
Info    : ImproveMesh 
Info    : Total badness = 1264.88 
Info    : Total badness = 1258.71 
Info    : SwapImprove  
Info    : 50 swaps performed 
Info    : SwapImprove2  
Info    : 0 swaps performed 
Info    : ImproveMesh 
Info    : Total badness = 1172.43 
Info    : Total badness = 1119.49 
Info    : CombineImprove 
Info    : 3 elements combined 
Info    : ImproveMesh 
Info    : Total badness = 1092.45 
Info    : Total badness = 1087.72 
Info    : SplitImprove 
Info    : badmax = 20.4893 
Info    : 0 splits performed 
Info    : ImproveMesh 
Info    : Total badness = 1087.72 
Info    : Total badness = 1087.73 
Info    : SwapImprove  
Info    : 22 swaps performed 
Info    : SwapImprove2  
Info    : 0 swaps performed 
Info    : ImproveMesh 
Info    : Total badness = 1053.89 
Info    : Total badness = 1029.86 
Info    : CombineImprove 
2024-04-11 13:36:06,063 [927] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.
Info    : 2 elements combined 
Info    : ImproveMesh 
Info    : Total badness = 1011.5 
Info    : Total badness = 1010.13 
Info    : SplitImprove 
Info    : badmax = 13.7339 
Info    : 0 splits performed 
Info    : ImproveMesh 
Info    : Total badness = 1010.13 
Info    : Total badness = 1010.07 
Info    : SwapImprove  
Info    : 13 swaps performed 
Info    : SwapImprove2  
Info    : 0 swaps performed 
Info    : ImproveMesh 
Info    : Total badness = 997.021 
Info    : Total badness = 988.982 
Info    : Optimizing volume 4
Info    : CalcLocalH: 187 Points 469 Elements 370 Surface Elements 
Info    : Remove Illegal Elements 
Info    : 169 illegal tets 
Info    : SplitImprove 
Info    : badmax = 13.4784 
Info    : 19 splits performed 
Info    : SwapImprove  
Info    : 20 swaps performed 
Info    : SwapImprove2  
Info    : 1 swaps performed 
Info    : 115 illegal tets 
Info    : SplitImprove 
Info    : badmax = 97.9023 
Info    : 15 splits performed 
Info    : SwapImprove  
Info    : 7 swaps performed 
Info    : SwapImprove2  
Info    : 0 swaps performed 
Info    : 82 illegal tets 
Info    : SplitImprove 
Info    : badmax = 97.9023 
Info    : 15 splits performed 
Info    : SwapImprove  
Info    : 1 swaps performed 
Info    : SwapImprove2  
Info    : 0 swaps performed 
Info    : 60 illegal tets 
Info    : SplitImprove 
Info    : badmax = 101.081 
Info    : 14 splits performed 
Info    : SwapImprove  
Info    : 1 swaps performed 
Info    : SwapImprove2  
Info    : 1 swaps performed 
Info    : 32 illegal tets 
Info    : SplitImprove 
Info    : badmax = 193.337 
Info    : 10 splits performed 
Info    : SwapImprove  
Info    : 1 swaps performed 
Info    : SwapImprove2  
Info    : 0 swaps performed 
Info    : 4 illegal tets 
Info    : SplitImprove 
Info    : badmax = 193.337 
Info    : 1 splits performed 
Info    : SwapImprove  
Info    : 0 swaps performed 
Info    : SwapImprove2  
Info    : 1 swaps performed 
Info    : 0 illegal tets 
Info    : Volume Optimization 
Info    : CombineImprove 
Info    : 17 elements combined 
Info    : ImproveMesh 
Info    : Total badness = 1473.71 
Info    : Total badness = 1375.72 
Info    : SplitImprove 
Info    : badmax = 51.5817 
Info    : 1 splits performed 
Info    : ImproveMesh 
Info    : Total badness = 1379.26 
Info    : Total badness = 1367.26 
Info    : SwapImprove  
Info    : 61 swaps performed 
Info    : SwapImprove2  
Info    : 0 swaps performed 
Info    : ImproveMesh 
Info    : Total badness = 1274.47 
Info    : Total badness = 1215.68 
Info    : CombineImprove 
Info    : 8 elements combined 
Info    : ImproveMesh 
Info    : Total badness = 1133.58 
Info    : Total badness = 1125.13 
Info    : SplitImprove 
Info    : badmax = 21.9556 
Info    : 1 splits performed 
Info    : ImproveMesh 
Info    : Total badness = 1130.48 
Info    : Total badness = 1129.19 
Info    : SwapImprove  
Info    : 21 swaps performed 
Info    : SwapImprove2  
Info    : 0 swaps performed 
Info    : ImproveMesh 
Info    : Total badness = 1101.38 
Info    : Total badness = 1081.27 
Info    : CombineImprove 
Info    : 8 elements combined 
Info    : ImproveMesh 
Info    : Total badness = 1008 
Info    : Total badness = 1003.95 
Info    : SplitImprove 
Info    : badmax = 12.3103 
Info    : 0 splits performed 
Info    : ImproveMesh 
Info    : Total badness = 1003.95 
Info    : Total badness = 1003.69 
Info    : SwapImprove  
Info    : 13 swaps performed 
Info    : SwapImprove2  
Info    : 0 swaps performed 
Info    : ImproveMesh 
Info    : Total badness = 996.234 
Info    : Total badness = 989.321 
Info    : Done optimizing mesh (Wall 0.19447s, CPU 0.197986s)
Info    : 772 nodes 4264 elements
Info    : Writing 'geometry/lv_ellipsoid.msh'...
Info    : Done writing 'geometry/lv_ellipsoid.msh'
2024-04-11 13:36:06,184 [927] DEBUG    UFL_LEGACY: Blocks of each mode: 
2024-04-11 13:36:07,024 [927] DEBUG    UFL_LEGACY: Blocks of each mode: 
2024-04-11 13:36:07,679 [927] WARNING  cardiac_geometries.fibers._utils: Unable to write XDMF file: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to find element family.
*** Reason:  Element Quadrature not yet supported.
*** Where:   This error was encountered inside XDMFFile.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  fd662efc1c0ddefa341c1dac753f717f6b6292a8
*** -------------------------------------------------------------------------
# Create the material
material_parameters = pulse.Guccione.default_parameters()
material_parameters["C"] = 10.0
material_parameters["bf"] = 1.0
material_parameters["bfs"] = 1.0
material_parameters["bt"] = 1.0
material = pulse.Guccione(parameters=material_parameters)
# Define Dirichlet boundary. Fix the base_spring
def dirichlet_bc(W):
    V = W if W.sub(0).num_sub_spaces() == 0 else W.sub(0)
    return DirichletBC(
        V,
        Constant((0.0, 0.0, 0.0)),
        geometry.ffun,
        geometry.markers["BASE"][0],
    )
# Traction at the bottom of the beam
lvp = Constant(0.0)
neumann_bc = pulse.NeumannBC(traction=lvp, marker=geometry.markers["ENDO"][0])
# Collect Boundary Conditions
bcs = pulse.BoundaryConditions(dirichlet=(dirichlet_bc,), neumann=(neumann_bc,))
# Create problem
problem = pulse.MechanicsProblem(geometry, material, bcs)
# Solve problem
pulse.iterate.iterate(problem, lvp, 10.0, initial_number_of_steps=200)
2024-04-11 13:36:08,064 [927] INFO     pulse.iterate: Iterating to target control (f_166)...
2024-04-11 13:36:08,065 [927] INFO     pulse.iterate: Current control: f_166 = 0.000
2024-04-11 13:36:08,065 [927] INFO     pulse.iterate: Target: 10.000
2024-04-11 13:36:08,150 [927] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-04-11 13:36:08,151 [927] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-04-11 13:36:08,151 [927] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-04-11 13:36:08,152 [927] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-04-11 13:36:08,856 [927] DEBUG    UFL_LEGACY: Blocks of each mode: 
  99	full
2024-04-11 13:36:09,063 [927] DEBUG    UFL_LEGACY: Blocks of each mode: 
  18	full
2024-04-11 13:37:06,773 [927] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-04-11 13:37:06,774 [927] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-04-11 13:37:06,774 [927] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-04-11 13:37:06,775 [927] INFO     UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-04-11 13:37:06,981 [927] DEBUG    UFL_LEGACY: Blocks of each mode: 
  10	full
2024-04-11 13:37:07,170 [927] DEBUG    UFL_LEGACY: Blocks of each mode: 
  3	full
([Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 115), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 182),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 115), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 212),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 115), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 237),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 115), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 262),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 115), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 287),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 115), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 312),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 115), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 337),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 115), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 370),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 115), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 411),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 115), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 460),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 115), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 509),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 115), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 550),
  Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 115), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 583)],
 [Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 180),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 210),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 235),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 260),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 285),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 310),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 335),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 368),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 409),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 458),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 507),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 548),
  Coefficient(FunctionSpace(None, FiniteElement('Real', None, 0)), 581)])
# Get displacement and hydrostatic pressure
u, p = problem.state.split(deepcopy=True)
endo_apex_marker = geometry.markers["ENDOPT"][0]
endo_apex_idx = geometry.vfun.array().tolist().index(endo_apex_marker)
endo_apex = geometry.mesh.coordinates()[endo_apex_idx, :]
endo_apex_pos = endo_apex + u(endo_apex)
print(
    ("\n\nGet longitudinal position of endocardial apex: {:.4f} mm" "").format(
        endo_apex_pos[0],
    ),
)
Get longitudinal position of endocardial apex: -26.5405 mm
epi_apex_marker = geometry.markers["EPIPT"][0]
epi_apex_idx = geometry.vfun.array().tolist().index(epi_apex_marker)
epi_apex = geometry.mesh.coordinates()[epi_apex_idx, :]
epi_apex_pos = epi_apex + u(epi_apex)
print(
    ("\n\nGet longitudinal position of epicardial apex: {:.4f} mm" "").format(
        epi_apex_pos[0],
    ),
)
Get longitudinal position of epicardial apex: -28.1947 mm
V = dolfin.VectorFunctionSpace(geometry.mesh, "CG", 1)
u_int = interpolate(u, V)
mesh = Mesh(geometry.mesh)
dolfin.ALE.move(mesh, u_int)
fig = plot(geometry.mesh, color="red", show=False)
fig.add_plot(plot(mesh, opacity=0.3, show=False))
fig.show()